package fem2.material.cohesive_law;

import fem2.Element;
import fem2.Variable;
import fem2.enu.VariableName;
import fem2.jlinalg.DiagonalMatrix;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import math2.MathUtilities;

/**
 * Mode I cohesive law based on Sancho paper: An embedded crack model for finite
 * element analysis of concrete fracture.<br>
 * Note: this class works with constant opening only
 * 
 * @author hbui
 * 
 */
public abstract class Mode1SanchoCohesiveLaw extends CohesiveLaw {

	// private Map<Element, Boolean> isLoading = new TreeMap<Element,
	// Boolean>();

	public Mode1SanchoCohesiveLaw(double Gf, double ft) {
		super(Gf, ft);
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLaw(double x);

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLawDerivative(double x);

	@Override
	public IVector computeJumpTraction(Element e, IVector w) {
		int dim = e.getMeshPart().getDimension();
		IVector t;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		if (wtilde < MathUtilities.tol6) {
			double ft = computeSofteningLaw(0);
			t = new ArrayVector(dim);
			t.set(0, ft);
			return t;
		}

		/*
		 * compute traction by checking for loading / unloading
		 */
		t = (IVector) w.clone();
		double wt_old = getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP).getValue();

		if (wtilde < wt_old) {

			/*
			 * unloading
			 */
			// double wt_old =
			// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP,
			// -1);
			double coeff = computeSofteningLaw(wt_old) / wt_old;
			BLAM.scale(coeff, t);

			// Debugger.watch(e.getId() + ", wt_old (t) = " + wt_old);

		} else {

			/*
			 * loading
			 */
			double coeff = computeSofteningLaw(wtilde) / wtilde;
			BLAM.scale(coeff, t);

			// Debugger.watch(e.getId() + ", wtilde (t) = " + wtilde);

			getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP).setValue(wtilde);

		}

		return t;
	}

	@Override
	public IMatrix computeJumpStiffness(Element e, IVector w) {
		int dim = e.getMeshPart().getDimension();
		IMatrix C;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		if (wtilde < MathUtilities.tol12) {
			C = new DiagonalMatrix(dim);
			double df = computeSofteningLawDerivative(0);
			C.set(0, 0, df);
			C.set(1, 1, 1e10);
			if (dim == 3) {
				C.set(2, 2, 1e10);
			}
			return C;
		}

		/*
		 * checking for loading / unloading
		 */
		double wt_old = getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP).getValue();
		if (wtilde < wt_old) {
			/*
			 * unloading
			 */
			C = new DiagonalMatrix(dim);
			// double wt_old =
			// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP,
			// -1);
			double coeff = computeSofteningLaw(wt_old) / wt_old;
			// C.set(0, 0, coeff);
			// C.set(1, 1, coeff);
			DiagonalMatrix.set(C, coeff);

			// Debugger.watch(e.getId() + ", wt_old (t) = " + wt_old);

		} else {
			/*
			 * loading
			 */
			C = new Array2DMatrix(dim, dim);
			double f = computeSofteningLaw(wtilde);
			double df = computeSofteningLawDerivative(wtilde);
			if (dim == 2) {
				double wn = w.get(0);
				double wm = w.get(1);
				double tmp1 = f / wtilde;
				double tmp2 = Math.pow(wn / wtilde, 2);
				double tmp3 = Math.pow(wm / wtilde, 2);
				double tmp4 = wn * wm / (wtilde * wtilde);
				double c00 = tmp1 * tmp3 + df * tmp2;
				double c01 = -tmp1 * tmp4 + df * tmp4;
				double c11 = tmp1 * tmp2 + df * tmp3;
				C.set(0, 0, c00);
				C.set(0, 1, c01);
				C.set(1, 0, c01);
				C.set(1, 1, c11);
			} else if (dim == 3) {
				double wn = w.get(0);
				double wm1 = w.get(1);
				double wm2 = w.get(2);
				double tmp1 = f / wtilde;
				double tmp2 = Math.pow(wn / wtilde, 2);
				double tmp3 = Math.pow(wm1 / wtilde, 2);
				double tmp4 = Math.pow(wm2 / wtilde, 2);
				double tmp5 = wn * wm1 / (wtilde * wtilde);
				double tmp6 = wn * wm2 / (wtilde * wtilde);
				double tmp7 = wm1 * wm2 / (wtilde * wtilde);
				double c00 = tmp1 * (tmp3 + tmp4) + df * tmp2;
				double c01 = -tmp1 * tmp5 + df * tmp5;
				double c02 = -tmp1 * tmp6 + df * tmp6;
				double c11 = tmp1 * (tmp2 + tmp4) + df * tmp3;
				double c12 = -tmp1 * tmp7 + df * tmp7;
				double c22 = tmp1 * (tmp2 + tmp3) + df * tmp4;
				C.set(0, 0, c00);
				C.set(0, 1, c01);
				C.set(0, 2, c02);
				C.set(1, 0, c01);
				C.set(1, 1, c11);
				C.set(1, 2, c12);
				C.set(2, 0, c02);
				C.set(2, 1, c12);
				C.set(2, 2, c22);
			}

			// Debugger.watch(e.getId() + ", wtilde (C) = " + wtilde);
			getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP).setValue(wtilde);
		}

		return C;
	}

	@Override
	public void Initialize(Element e) {
		Variable v = new Variable(0);
		v.registerScalarVariable(VariableName.EQUIVALENT_JUMP, false);
		setVariable(e, v);
		// isLoading.put(e, false);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		getVariable(e).save();
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
